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Abstract 

Skyrmions are spiral structures observed in thin films of certain mag- 
netic materials [T]. Of the phases allowed by the crystalline symmetries 
of these materials [2] only the hexagonally packed phases (SCh) has been 
observed. Here the melting of the SCh phase is investigated using Monte 
Carlo simulations. In addition to the usual measure of Skyrmion density 
chiral charge, a morphological measure is considered. In doing so it is 
shown that the low temperature reduction in chiral charge is associated 
with a change in Skyrmion profiles rather than Skyrmion destruction. At 
higher temperatures the loss of six fold symmetry is associated with the 
appearance of elongated Skyrmions that disrupt the hexagonal packing. 
Pacs :12.39.Dc, 75.70.-i 



1 1. School of Physics, The University of Western Australia, 35 Stirling Hwy, Crawley 6009 
Australia 

ambrom01@student.uwa.edu.au 

2 2. SUPA School of Physics and Astronomy, University of Glasgow, Glasgow G12 8QQ 
United Kingdom 
Robert . Stamps@glasgow. ac.uk 



1 Introduction 

Following early work by Pauli a variety of non-linear sigma models were sug- 
gested as models for baryons [21 2] [5]. One of the first soluble models was 
found by Skyrme who considered rotational variables [6] and the resultant unit 
field spiral solitons are referred to as Skyrmions. Two dimensional Skyrmions 
appeared as the fundamental excitation from a spin polarized two dimensional 
electron gas. It was suggested that defects could localize Skyrmions giving rise to 
the fractional quantum hall effect [7J [H [5] • Observations in GaAs quantum well 
systems using NMR 1Q\ and magnetoabsorption spectroscopy [11] confirmed the 
presence of quasi-particles with charges consistent with theoretical predictions. 
More recently measurements have been made in GaAs using NMR relaxation 
|12j . spin wave absorption |13) and microwave absorption [14] . that are con- 
sistent with the presence of a predicted Skyrmion lattice |15j . Skyrmions were 
first proposed as a possible magnetic spin texture by Bogdanov and Hubert who 
showed that Dzyaloshinskii-Morya exchange terms could lead to stabilization of 
Skyrmion crystals in chiral magnets |16j . Recently the real space measurements 
of the two dimension analogue these spiral structures has been made at low 
temperature in thin films of Feo.5Coo.5Si [Tl 117) and close to room temperature 
in FeGc[18j, along with the measurement of spin torque effects at very low cur- 
rent densities |19j . these measurements have ignited interest in Skyrmions as a 
possible candidate for magnetic storage [2D]. 




In two dimensions Skyrmions exist as a vortex structure modulated by a chang- 
ing perpendicular component. Consider a field of unit length s(x) with x two 
dimensional and s(x) taking the usual polar representation [^] Away from an 
isolated Skyrmion the field is described by zenith angle 9 = tt. Taking radial 
coordinates (p, if)) for x a Skyrmion can be described as a vortex in the az- 
imuthal angle if>(p,ip) — ip — ir/2 modulated by a radial varying zenith angle 
9(p,ijj) = 9(p) such that 9(0) — and 9(R) = tt. An example of close packed 
Skyrmions projected into the plane is shown in figure [I] (c) where gray(black) 
arrows indicate spins with perpendicular component are pointing up (down). In 
order for such an excitation to be stable one requires a Hamiltonian with cither 
fourth order derivatives of the field [H] or terms lacking inversion symmetry [TB]. 
Lack of inversion symmetry can be found in a variety of crystal structures, in 
particular B20 crystal compounds with crystallographic point group 23 in 
which Dzyaloshinskii-Morya (DM) coupling is present. Yi et a.l [2] calculated 
the T = phase diagram as a function of external magnetic field and anisotropy 
for such a two dimensional magnet using Monte-Carlo simulation. Five different 
states were identified: The high field saturated ferromagnetic phase, a helical 
state (figure [IJa)) and three chiral states shown. The chiral crystal states were 
labeled according to their rotational symmetries: SC\ (figure [lib)) has a two 
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Figure 1: Projection of spins into the x-y plane, black arrows indicate that 
sf > while gray arrows indicate sf < A: Helix State; B: SC% Skyrmion; C: 
SC2 Skyrmion; D: SCh Skyrmion. 



SCh (figure [T] (c)) a six fold symmetry Q In what follows we restrict our atten- 
tion the the SCh type skyrmions. 

Analytically the ground state properties of this two dimensional case have been 
studied using a Landau-Ginzburg formalism |21j . The nature of the phase tran- 
sition has already been investigated using combinatorics [22]. Here the au- 
thors considered the striped helical state structure shown in figure [l] (a) as the 
ground state. In the presence of an external field perpendicular to the plane 
of the sample, one direction of perpendicular magnetization is favored. Since 
the stripe period is fixed by the ratio of Dzyaloshinskii-Morya and ferromag- 
netic exchange coupling, spins reverse inside stripes anti-aligned with the field, 
breaking remaining anti-aligned areas into finite length stripes. Unlike the one 
dimensional periodicity of the ground state (giving zero chiral density) the ends 
of the terminating stripes form a 'meron'; a chiral half spiral. The Skyrmion 
crystal state SCh show in figure [ljc) is the state with maximum possible meron 
density. By considering the merons as particles with chiral charge, a free energy 
can be written and a qualitatively correct B-T phase diagram produced. While 
the assumptions of this model are reasonable low temperature, the possibility 
of Skyrmion profiles changing with temperature cannot be described. One ex- 
pects, that as temperature is increased, thermal fluctuations will dominate and 

5 This notation should not be confused with the notation some authors in which Skyrmion 
states are labeled as either SkX or SkG refering to whether the SCh nas formed a close packed 
crystal (SkX) or are freely moving (SkG) 
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the particle description will no longer be valid. 

Current Monte Carlo simulations have focused on obtaining phase diagrams, 
with phases identified according to symmetry properties and chiral density 
[171 [2]. As a system with textured phase, close packed Skyrmions represent 
three types of order: a net chiral charge, a six fold symmetry and a net magne- 
tization. Here we use Monte Carlo simulations to investigate the melting of the 
SCh phase, focusing on the mechanism through which these various orders are 
destroyed. 

2 Theory and Method 

Bak and Jensen proposed a continuum Hamiltonian [23 to explain observations 
of helical structure in MnSi [2HE3I26] and FeGe [27]. The model was adapted 
to a discrete lattice by Yi et al. [2] and we begin with their Hamiltonian. In 
two dimensions and in our notation, it reads 

E = - J/2 ^ &i ' 9j ~ K / 2 X (Si+ax - Si-az)) ■ X + (siX (s i+ay - Si-aj))) ' $) 

i i i 

(i) 

For brevity we use single subscripts to represent locations on a two dimensional 
square lattice with lattice constant a. The Sj are dimensionless spins at the 
vertices of this lattice representing the average magnetization of a small volume 
Vj = a x a x t centered on vertex i, with t this thickness of the film. Si +ax 
and Si-ax represent the closest spins to in the x direction (with analogous 
notation in the y direction) and Xwjj) indicates a sum over all pairs of nearest 
neighbors. K is the strength of the Dzyaloshinskii-Morya exchange coupling, 
the sign of which determines the handedness of the Skyrmions; J is the strength 
of the isotropic exchange coupling; H is the strength of the applied magnetic 
field and A\ and A 2 are anisotropies. The ratio of exchange to Dzyaloshinskii- 
Morya strength determines the periodicity of Skyrmion lattice as P = da where 
dis[l[T7]: 



d = 2?r 



(2) 



Since the computation time to complete a Monte Carlo step scales with the 



Parameter 


Value 


K/J 

Ax/J 

A 2 /J 


2.45 
0.5 
0.5 



Table 1: Parameters used in the Monte Carlo simulation. 
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number of spins, simulating large systems is typically associated with long com- 
putation times. In order to decrease this computation time GPU parallel pro- 
gramming is employed. The Hamiltonian includes only short range interactions 
so a parallel checkerboard type algorithm similar to those described by Weigel 
and Yavorskii [2"81 |2"9"] is used. 

2.1 Ground State Packing 

In the interest of comparison with previous results [T71 15] we choose the ma- 
terial parameters so that the expected dominant wave length will be d = 6. 
Parameters used are given in table [l] where dimensionless quantities are given 
by normalizing against J. In what follow temperature and energy will be given 
as T = {ksT^/J and energy as £ = E/J. The magnetic field is fixed at 
H/J = 1.875 to ensure that at T = the system forms close packed Skyrmions. 
Some care should be taken when choosing the finite dimension (L) of the spin 
lattice for these calculations. For Skyrmions the ground state is periodic, con- 
sisting of hexagonally close packed Skyrmions, hence the repeated unit cell of 
the crystal lattice is not square. In order to tile a pattern with six fold symme- 
try one usually selects a unit cell that is a parallelogram defined by the vectors 
vi = (d,0) and i>2 = (d/2, v3rf/2), where d is the spacing between six fold 
symmetry axes (in this case Skyrmion cores). The resulting parallelogram is 
indicated in figure [2] This parallelogram unit cell is space filling and can be 




Figure 2: Left: Unit cells that span a plane with hexagonal symmetry; Right: 
the dimensions of a rectangular unit cell that can be used to tile a finite system 
with periodic boundary conditions. 

used to cover an infinite plane. 

Our simulation uses a square lattice spanning a finite region with periodic 
boundary conditions. In order to tile a space with periodicity in the x and 
y directions we select a rectangle to act as the unit cell, also shown in figure 
[2l This unit cell is then a rectangle with width d\ = d and height e?2 = V3d. 
The above parameters for which d — 6 correspond to a unit cell of dimensions 
6 x 10.3923. If one picks a system size with periodic boundary conditions and 
length as multiples of d, the resultant state will be influenced in one direction 
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by the forced periodicity. When deciding the system size one should attempt 
to find a system size that is a multiple of d\ and d 2 . With d 2 irrational this 
is not possible. Instead we consider the two lengths that define the unit cell 
di = 6 and d 2 — 6\/3- In order to measure the amount of mismatch between 
the system size L and the unit cell we define a function: 



M(L, di,d 2 ) =min(L - d 2 [L/d 2 \ ,d 2 - L + d 2 [L/d 2 \ ) 
+min(L - d t [L/d x \ ,d\ — L + d\ \L/d\\ ) 



(3) 



where each term in equation (J3| is a measure of how far L is from an exact 
multiple of d\ or d 2 . In the case of the hexagonal Skyrmion lattice where d 2 is 
irrational and L is restricted to be a integer, the minima are not known a priori. 
In order to select a reasonable size, M(L, 6, 6-\/3) was calculated for L G [1, 100]. 
As shown in figure [3j the minimum value is found to be L = 42. 




40 60 
System Size (L) 



100 



Figure 3: Log plot of M(L, 6, 6^/3) showing the minimum at L = 42. This 
minimum corresponds corresponds to the tiling that best minimizes finite size 
effects introduced by periodic boundary conditions. 



2.2 Analysis 

The close packed Skyrmion ground state consists of two types of order, a chiral 
charge associated with the creation of spiral structures and a long range order 
associated with the hexagonal packing. In addition to the Skyrmion crystal 
order the system has a net magnetization in response to the applied magnetic 
field. 

Two techniques are employed in order to measure the number of Skyrmions in 
a state: The first is to consider the topological charge Q of a Skyrmion , which 
is the sum of the local chiral charge over the area of a single Skyrmion [22] . The 
charge density is: 

Pi ~j Si ' ((^-fax ax) X (s^+ay ^i— ay)) (4) 
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In order to define the charge Q of a Skyrmion, one calculates the total charge 
density of the ground state and divides by the number of Skyrmions. In order 
to calculate the Skyrmion number at some finite temperature one divides the 
total charge by Q. 

This method of counting Skyrmions assumes that the Skyrmion charge remains 
constant (i.e that the Skyrmion profile is not temperature dependent). For 
comparison, we define a second measure of Skyrmion number. Since the core of 
a Skyrmion points against the applied magnetic field one can binarize a state by 
applying 0(l/2(s^.l + 1) — T r ) (0 is the Heaviside Theta function), where T r is 
some threshold between zero and one. In doing so one identifies spins with s^.z 
anti-parallel to the magnetic field as Skyrmion cores. One can then calculate 
the connected components of the resultant state. To do this the eight neighbors 
and next nearest neighbors to a given spin are considered connected if they are 
equal. Counting the total number of connected components gives a measure of 
the Skyrmion number. By comparing these two measures of Skyrmion number 
one can distinguish between reduction of topological charge due to Skyrmion 
destruction and alteration of Skyrmion profile. 

In order to examine the long range order of a state the Fourier transform Sk = 
Y]j Sj exp(«fc • j) is taken. When the the system is close packed |<Sfc| will have six 
satellite peaks corresponding to the hexagonal symmetry of the ground state. 

3 Results 

For each system size to be investigated the ground state is first calculated by 
starting the system in a random configuration and then reducing the tempera- 
ture to zero over 10 7 MC steps. The ground state for L = 42 in shown in figure 
[4] along with the intensity profile of the Fourier transform | Sk \ • By selecting the 
system size L to match the dimensions of the rectangular unit cell described 
above, it is possible to fit an integer number of unit cells into the square array. 
For finite temperature results the system is evaluated at a constant temperature 
starting from the ground state with the first 10 7 MC steps disregarded to allow 
the system to equilibrate. An ensemble of 100 states is calculated and 280 MC 
steps are taken between subsequent states to ensure that each state is selected 
independently. In figure [5] the order parameters of the system are shown indi- 
cating the low temperature transition into the Skyrmion state. The persistence 
of the magnetic response at temperatures greater than those required for the 
destruction of chiral order is seen through the magnetization M z and the cone 
angle. The loss of chiral structure occurs around T = 3, but the elevated cone 
angle and magnetic order persist at approximately T = 8. The energy variance 
is double peaked indicating two broad transitions. At high temperatures the 
entropy dominates and there is no order. Below T = 10 there exists a region in 
which the cone angle increases from its high temperature value giving rise to a 
non zero magnetization. The external field ensures that as the temperature is 
decreased the magnetization increases. Below approximately T = 2.5 there is 
a transition into the Skyrmion state characterized by the sharp increase in the 
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Figure 4: The ground state spin configuration (top) and the intensity of the 
Fourier transform (bottom) showing peaks characteristic of the six fold transla- 
tional symmetry of the hexagonal lattice. 
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Figure 5: Clockwise from top left: The perpendicular magnetization, the chiral 
charge, the cone angle and the energy variance. 



strength chiral density (negative values are indicative of left handed structure). 
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Since the creation of Skyrmions prevents the system reaching saturation we note 
that the magnetization remains well below the saturation value. 



3.1 Loss of six fold order 

In order to examine the loss of long range order we focus on the low tempera- 
ture regime. In figure[6]the Skyrmion number is calculated using the total chiral 
charge and the method of counting cores described previously. While the chiral 
measure is reduced at all finite temperatures the number of reversed regions 
remains constant below about T = 0.3. This suggests that the initial loss of 
chiral charge is due to thermal distortion of the Skyrmion profiles rather than 
destruction of the Skyrmions themselves. In figure [7] an example state from 
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Figure 6: Two measures of Skyrmion packing as a function of temperature: the 
connected reversed components (gray circles) and the ratio of total chiral charge 
to the charge of one Skyrmion in the ground state p/Q (black circles). 

ensembles at different temperatures is show along with \Sk\, between T = 0.1 
and T — 2.5. At T = 0.4 the strong six fold symmetry is compromised with 
peaks beginning to smear together. At T = 0.5 the peaks have smeared into a 
circle, indicating that a while a dominant length scale still exists sixfold transla- 
tional symmetry is lost. At this temperature one can simultaneously observe the 
appearance of elongated regions of reversed magnetization. As temperature is 
further increased any order in the perpendicular components of s*j is destroyed. 
An enlarged example of a state at T = 0.5 is shown in figure [8] In addition to 
the perpendicular magnetic order, in plane magnetic order is shown as arrows. 
In addition to disrupting the close packing of Skyrmions these extended struc- 



To illustrate this, 
8] is replotted with 



tures also reduce the chiral density described in equation [4 
in figure [9] the longest of extended regions shown in figure 
information about perpendicular order omitted. Color is used to indicate the 
areas of highest chiral density, blue represents areas of low density and and red 
indicates areas of high density. Here only the ends of the structure contribute 
significantly. 
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Figure 7: The perpendicular components of si and the Fourier intensity | jS'fc | 
near the phase transition. Since the central (k = 0) Fourier peak will signifi- 
cantly larger than the satellite peaks at finite temperature, it is removed to give 
contrast. The Fourier plots then have their colors scaled with green representing 
zeros and and red representing the maximum value. Perpendicular components 
of Si are colored according to the legend above. Left column from top: T = 0.1, 
T = 0.2, T = 0.35, T = 0.4 Right column from top: T = 0.5, T = 0.8, T = 1.2, 
T = 2.5. 
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Figure 8: An example of the elongated structures associated with the destruction 
of long range order at T = 0.5. Here the perpendicular components of the spins 
are represented by color as in figure [7J in addition for spins parallel to the plane 
their direction is indicated with arrows. 




Figure 9: An example of the elongated structures associated with the destruction 
of long range order at T = 0.5. Spins parallel to the plane have their direction 
indicated with arrows. Areas of high chiral density are indicated in red. 

The presence of the low temperature (below T = 0.3) distortion of Skyrmion 
profile and the appearance elongated regions might indicate a change of the 
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L=39 






Figure 10: (\S k \) as function of size for T = 0.2(Top) and T = 0.4(Bottom). 
The plots are colored as in figure [7| 

dominant length scale. While the difficulty associated with finite size effects 
was minimized by considering the dimensions of the ground state unit cells, 
this doesn't necessarily ensure that the effects are negligible at finite tempera- 
tures. Several other comparable system sizes were investigated for T = 0.2 and 
T = 0.4. The resulting energies are given in table [2] where it is seen that at 
finite temperature the change in energy due to finite size effects is not signifi- 
cant. 



L 


<£>(T = 0.2) 


(£ >(T =0.4) 


39 


-2.31898 


-2.14408 


42 


-2.32090 


-2.14087 


45 


-2.32448 


-2.13814 


48 


-2.31905 


-2.13881 



Table 2: Ensemble energies for (£) at T = 0.2 and T = 0.4. 
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If the spacing 



(|S/c|) was also calculated and the results are shown in figure 
between Skyrmions was to change with temperature, one expects a small change 
in system size could stabilize six fold ordering at higher T. Systems with size L 
equal to 45 or 48 show the same general trend in which the six order is destroyed 
leaving peaks in the [1, 1] and [—1,1] directions. For L = 39 the system doesn't 
form six fold packing. These results are consistent with the above interpretation 
of the phase transition as a result of changing Skyrmion profile. 



4 Conclusions and Comments 

We have shown that since the chiral density changes smoothly over a large tem- 
perature range, it does not capture all of the information about the melting 
process from the SCh phase in chiral magnets. The Fourier intensity and the 
reversed connected components of the system reveal a low temperature reduc- 
tion in the chiral charge due to thermal distortion of the Skyrmion profiles. 
At higher temperatures there is a sharp loss of six fold order associated with 
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creation of elongated structures that disrupt ordering. Only the ends of these 
structures contribute to the chiral charge, further reducing the chiral density. 
In addition it has been shown that the nature of hexagonal close packing means 
that while any choice of simulation size will introduce finite size effects, judi- 
cious choice of system size can help alleviate these effects. 

Here the authors have limited the scope of investigation to a single choice 
of applied field, however, the techniques presented might offer insight for a 
wide variety of fields. Specifically there is the possibility that at higher fields 
the relative strength J2i H.Si might suppress the creation of extended regions 
of reversed magnetization. There exists also the possibility of analyzing the 
zero field ground state as stripes of alternating connected regions of perpen- 
dicular spins pointing in the positive and negative z direction, represented as 
0(l/2(±Si.z + 1) — T r ). Similarly the SC\ and SC2 phase might be analyzed 
as a checkerboard alternating cores. 

Interestingly the system maintains a net chiral charge at temperatures at which 
any long range ordering is destroyed which could have consequences for the 
calculation of hall conductivity at these temperatures [2] . 
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